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Abstract 

We give a full description of the code NBODY2 for direct integration of the gravita- 
tional iV-body problem. The method of solution is based on the neighbour scheme 
of Ahmad & Cohen (1973) which speeds up the force calculation already for quite 
■ modest particle numbers. Derivations of all the relevant mathematical expressions 

are given, together with a detailed discussion of the algorithms. The code may be 
used to study a wide variety of self-consistent problems based on a small softening 
i ' of the interaction potential. 
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iV-body simulations provide a useful tool for exploring a wide range of as- 
tronomical problems. By now, direct summation methods have been used to 
study such diverse topics as planetary formation, few-body scattering, star 
cluster dynamics, violent relaxation, tidal disruption of dwarf galaxies, inter- 
acting galaxies, compact groups and clusters of galaxies, as well as cosmologi- 
cal modelling. The code nbody2 described here is suitable for relatively small 
iV-body systems (20 < N < 10 4 ). However, it does not contain any special 
treatments of close encounters and should therefore only be used for problems 
where the effect of hard binaries can be neglected (a variety of regularization 
methods are discussed in Aarseth 1985). 
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2 Integration Method 



In this Section, we provide all the relevant mathematical expressions and dis- 
cuss the basic principles of iV-body integration. 



2.1 Difference Formulation 



We write the equation of motion for a particle of index % in the form 



N I 

mi r, ; - r. 



G J^,(4+V /2 - (1) 



Here G is the gravitational constant and the summation is over all the other 
particles of mass irij and coordinates r,,-. The introduction of a softening pa- 
rameter e prevents a force singularity as the mutual separation — > 0, thereby 
making the numerical solutions well behaved. For convenience, we adopt scaled 
units in which G = 1 and define the left-hand side of (1) as the force per unit 
mass, F, omitting the subscript. The present difference formulation is based 
on the notation of Ahmad & Cohen (1973, hereafter AC) and follows closely 
an earlier treatment (Aarseth 1985). 

Given the values of F at four successive past epochs t 3 , t 2 , t±, t , with t the 
most recent, we write a fourth-order fitting polynomial at time t as 

F t = (((D 4 (t - t 3 ) + D 3 ) (t - t 2 ) + D 2 ) (t - ti) + D 1 )^ - t ) + F„. (2) 

Using compact notation, the first three divided differences are defined by 

V k Kt k ] = Dfc ' 1[t0 ' tfc : l] " Dfc ' 1[tl ' tfc] ! (k = 1, 2, 3) (3) 

to — tk 



where D° = F and square brackets refer to the appropriate time intervals 
(D 2 [t!,t 3 ] is evaluated at ti). The term D 4 is defined similarly by D 3 [t,t 2 ] 
and D 3 [t , £3]. Conversion of the force polynomial into a Taylor series provides 
simple expressions for integrating coordinates and velocities. Equating terms 
in the successive time derivatives of (2) with an equivalent Taylor series and 
setting t — t yields the force derivatives 



F« = ((D 4 t 3 + D 3 K 2 + D 2 ) t[ + D 1 

F( 2 ) = 2! {D'it'A + t' 2 t' 3 + f/3) + D 3 (t; + 4) + D 2 ) 
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F< 3 > = 3! (D 4 (t; + 4 + Q + D 3 ) 

F (4) = 4! D 4 , (4) 



where t' k = to — t/~. Equation (4) is mainly used to obtain the Taylor series 
derivatives at t — to, when the fourth difference is not yet known. Thus the 
contribution from D 4 to each order is only added at the end of an integration 
step. This so-called 'semi-iteration' gives increased accuracy at little extra cost 
(on scalar machines) and no extra memory. 

We now describe the initialization procedure, assuming one force polynomial. 
From the initial conditions rrij, rj, Vj, the respective Taylor series derivatives 
are formed by successive differentiations of (1). Introducing the relative co- 
ordinates R = Yi — Yj and the relative velocities V = — Vj, all pair- wise 
interaction terms in F and F^ are first obtained by 



rrijR 

ij R* 

Fg> = -!^-3aF y , (5) 

with a = R • V/R 2 . The total contributions are obtained by summation over 
N. Next, the mutual second- and third-order terms are formed from 



F (2) = 

13 



3b F 



l 3 



raj (F, 



(i) 



R 3 



- 9aF ( ff - 96Fg } - 3cF ij? 



(6) 



with 



6=(|) a + R ' ( ^T Fj) +a a 

c = ^ + ~Jp — +a(3b-Aa 2 ). (7) 

A second double summation then gives the corresponding values of F^ 2 ^ and 
F^ 3 ) for all particles. This pair-wise boot-strapping procedure provides a con- 
venient starting algorithm, since the extra cost is usually small. 

Appropriate initial time-steps Atj are now determined, using the general cri- 
terion discussed in the next subsection. Setting to = 0, the backward times are 
initialized by t^ = —kAti (k = 1, 2, 3). Inversion of (4) to third order yields 
starting values for the divided differences, 
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D 1 



( i F (3) t ; _ ' F (2) )t[ + F (i) 



D 2 



- l -F^ (t[+t' 2 ) + V> 



D 3 



-F (3) . 
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(8) 



2.2 Individual Time-Steps 

Stellar systems are characterized by a range in density which gives rise to 
different time-scales for significant changes of the orbital parameters. In order 
to exploit this feature, and economize on the expensive force calculation, each 
particle is assigned its own time-step which is related to the orbital time- 
scale. Thus the aim is to ensure the convergence of the force polynomial (2) 
with the minimum number of force evaluations. Since all interactions must be 
added consistently in a direct integration method, it is necessary to include a 
temporary coordinate prediction of the other particles. However, the additional 
cost of low-order predictions still leads to a significant saving since this permits 
arbitrarily large time-step ratios. 

Following the polynomial initialization discussed above, the integration cycle 
itself begins by determining the next particle (i) to be advanced. Thus in 
general, % = min., (tj + At,), where tj is the time of the last force evaluation. 
It is convenient to define the present epoch (or 'global' time) by t = U + At,, 
rather than relating it to the previous value. 

The individual time-step scheme (Aarseth 1963) uses two types of coordinates 
for each particle. We define primary and secondary coordinates, r and r t , 
evaluated at t and t, respectively, where the latter are derived from the former 
by the predictor. In the present treatment where high precision is not normally 
required, we predict coordinates to order F^ 1 ) by 



where F^ = |FW, F = §F and bt'- =t-t j (with St'j < Atj). The coordinates 
and velocities of particle i are then improved to order F^ 3 ^ by standard Taylor 
series integration (cf. (4)), whereupon the current force may be obtained by 
direct summation. At this stage the four times tk are updated to be consistent 
with the definition that t denotes the time of the most recent force evalua- 
tion. New differences are now formed (cf. (3)), including D 4 . Together with 
the new F^ 4 \ these correction terms are combined to improve the current co- 
ordinates and velocities to highest order, whereupon the primary coordinates 
are initialized by r = r t . 



((F^St'j + F)5t' j + v)5t' j + r , 



(9) 



4 



New time-steps are assigned initially for all particles and at the end of each 
integration cycle for particle i. We adopt the composite criterion 

_ ( ifflffligf) ) 1 '' (10) 



where r/ is a dimensionless accuracy parameter. For this purpose, only the 
last two terms of the first and second force derivatives in (4) are included. 
This expression ensures that all the force derivatives play a role and is also 
well defined for special cases (i.e. starting from rest or |F| ~ 0). Although 
successive time-steps normally change smoothly, it is prudent to restrict the 
growth by a stability factor (presently 1.2). 

In summary, the scheme requires the following 30 variables for each particle: 
m, r , r t , v , F, F^, D 1 , D 2 , D 3 , At, t , t±, t 2 , t 3 . It is also useful to employ 
a secondary velocity, v t , for dual purposes. 



2.3 Neighbour Scheme 



Having introduced the basic tools for direct integration, we now turn to the 
AC neighbour scheme. Here the main idea is to reduce the effort of evaluating 
the force contribution from distant particles by employing two polynomials 
based on separate time-scales. Splitting the total force on a given particle into 
an irregular and a regular component, 

F F; rr -|- F re g , ( 1 1 J 



we can replace the full N summation in (1) by a sum over the n nearest 
particles together with a prediction of the distant contribution. This procedure 
can lead to a significant gain in efficiency, provided the respective time-scales 
are well separated and n « N. 

To implement the AC scheme, we form a list for each particle containing all 
members inside a sphere of radius R s . In addition, we include any particles 
within a surrounding shell of radius 2 1 / 3 R s satisfying R - V < 0.1i? 2 /ATj, 
where AT, denotes the regular time-step. This ensures that fast approaching 
particles are selected from the buffer zone. 

The size of the neighbour sphere is modified at the end of each regular time- 
step when a total force summation is carried out. We have adopted a neighbour 
criterion mainly based on distance for convenience. However, some modifica- 
tion would be desirable for very large mass ratios (say > 100). A selection 
criterion based on the local number density contrast has proved itself for a 
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variety of problems, but may need modification for interacting subsystems 
or massive particles. To sufficient approximation, the local number density 
contrast is given by 



where n\ is the current membership and Rh is the half-mass radius. In order 
to limit the range, we adopt a predicted membership 

n p = n max (0.04C) 1/2 , (13) 



subject to n p being within [0.2n max , 0.9n max ], with n max denoting the maxi- 
mum permitted value. The new neighbour sphere radius is then adjusted using 
the corresponding volume ratio, which gives 

?0 id fM 1/3 



^new = R oia | _J> j _ ^ 



An alternative and simpler strategy is to stabilize all neighbour memberships 
on the same constant value n p = Hq (AC, Makino & Hut 1988). A number of 
refinements are also included as follows: 

• In order to avoid a resonance oscillation in R s , (rip/rii) 1 ^ is used in (14) if 
the predicted membership lies between the old and new value 

• R s is modified by a radial velocity factor outside the core 

• The volume factor is only allowed to change by 25 %, subject to a time-step 
dependent cut-off if AT; < 0.01 T cr (T cr is the crossing time) 

• If ni < 3 and the neighbours are moving outwards, the standard sphere R s 
is increased by 10 %. 

The gain or loss of particles is recorded when comparing the old and new 
neighbour list, following the re-calculation of the regular force. Regular force 
differences are first evaluated, assuming there has been no change of neigh- 
bours. This gives rise to the provisional new regular force difference 

"pmcw ("Pold rpncw\ xpold 

T-vl _ rcg \ r irr r irr ) r reg /ir-i 



where F°'g denotes the old regular force, evaluated at time T , and the net 
change in the irregular force is contained in the middle brackets. In the subse- 
quent discussions, regular times and time-steps are distinguished using upper 
case characters. All current force components are obtained using the predicted 
coordinates (which must be saved), rather than the corrected values based 



6 



on the irregular D 4 term since otherwise (15) would contain a spurious force 
difference. The higher differences are formed in the standard way, whereupon 
the regular force corrector is applied (if desired). 

A complication arises because any change in neighbours requires appropriate 
corrections of both force polynomials; i.e. using the principle of successive 
differentiation of (11). The respective Taylor series derivatives (5) and (6) are 
accumulated to yield the net change. Each force polynomial is modified by first 
adding or subtracting the correction terms to the corresponding Taylor series 
derivatives (4) (without the D 4 term), followed by a conversion to standard 
differences using (8). 

Implementation of the AC scheme requires the following additional set of reg- 
ular variables: F reg , D 1 , D 2 , D 3 , AT, T , 7\, T 2 , T 3 , as well as the neighbour 
sphere radius R s and neighbour list L (size n max + 1). 



3 Program Structure 

This Section describes the structure of the code and defines the main input 
parameters and some of the options. 

3. 1 Routines 

The program contains some 2307 lines of Fortran statements and another 1203 
lines of comments and spaces. It is divided into 38 routines, many of which 
are optional and not required for standard simulations. All the routines are 
listed in Table 1, together with the main calling sequence and an outline of 
the purpose denotes an option). 

The program consists of four separate parts: input, output, initialization and 
integration. The basic routines are: start, input, data, scale (input); energy, 
core, output (output); nblist, fpolyi, fpoly2, steps (initialization); and 
intgrt, nbint, REGiNT, stepi (integration). The main calling sequence (start, 
output, intgrt) originates in routine main, whereas routine start acts as 
driver for both input and polynomial initialization and routine intgrt controls 
the subsequent program flow. 

To improve the layout and distinguish between code and text, we have adopted 
upper case characters for the former. At the 'microscopic' level, blank spaces 
are used to improve readability of the Fortran statements, where blocks of 
code inside Do-loops and logical /f's are indented. 
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Table 1 

Names of routines 



Routine 


Called by 


Description 


BLOCK 




Block data initialization 


BODIES 


OUTPUT 


Output of binaries and particles 6 k 9) 


CHECK 


OUTPUT 


Error control and restart 2 k 11) 


CMCORR 


OUTPUT 


Centre of mass corrections (# 18) 


COAL 


MAIN 


Inelastic two-body collision (# 12) 


CORE 


OUTPUT 


Core radius and density centre 8) 


CPUTIM 


OUTPUT k INTGRT 


CPU time in minutes (machine dependent) 


DATA 


START 


Generation of initial conditions 


DEFINE 


INPUT 


Definition of input parameters and counters 


ENERGY 


OUTPUT k SCALE 


Total energy 


ESCAPE 


OUTPUT 


Removal of escaping particles (# 13) 


FPCORR 


REGINT 


Derivative corrections of force polynomials 


FPOLYl 


START k COAL 


Force and first derivative 


FPOLY2 


START k COAL 


Second and third force derivatives 


INPUT 


START 


Main parameter input 


INSERT 


INTGRT 


Insert of particle index in sequential list 


INTGRT 


MAIN 


Decision-making and flow control 


LAGR 


OUTPUT 


Lagrangian radii and half-mass radius 7) 


MAIN 




Master control flow 


MODIFY 


MAIN 


Reading modified parameters at restart 


MYDUMP 


MAIN k INTGRT 


common save or copy (# 1 k 2) 


NBINT 


INTGRT 


Irregular integration 


NBLIST 


START k COAL 


Initialization of neighbour list 


OUTPUT 


MAIN 


Main output and data save 


RAN2 


DATA 


Portable random number generator 


REGINT 


INTGRT 


Regular integration 


REMOVE 


ESCAPE k COAL 


Removal of particle from common arrays 


SCALE 


START 


Scaling to new units 


SORT2 


INTGRT k LAGR 


Sequential sorting of array 


START 


MAIN 


Calling sequence for initial setup 


STEPI 


INTGRT 


Time-step for irregular force polynomial 


STEPS 


START k COAL 


Initialization of time-steps and differences 


SUBSYS 


START 


Initial subsystems for binary test (#17) 


VERIFY 


INPUT 


Validation of main input parameters 


XTRNL1 


REGINT k FPOLYl 


Force due to Plummer potential (# 15) 


XTRNL2 


REGINT k FPOLYl 


Force due to logarithmic potential (# 15) 


XVPRED 


OUTPUT k COAL 


Prediction of coordinates and velocities 


ZERO 


START 


Initialization of global scalars 



3.2 COMMON Variables 



In principle, the program can be used for any particle number, subject to 
available computer memory and CPU time. The size of all particle arrays is 
controlled by the parameter nmax, together with lmax which is used for the 
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neighbour lists. These parameters must be specified in the auxiliary header 
file params . h at compile time since most routines communicate via the global 
common blocks, rather than by general-purpose arguments. All the common 
variables are defined in a T^K file for convenience. 2 

The global labelled common blocks (header file common2.h) are organized in 
four separate parts denoted nbody, names, counts, params, which contain 
all the essential variables. This enables a calculation to be halted temporarily 
by saving the common blocks on disc or tape, whereupon a restart can be 
made directly from the high-order scheme. 

The convention of up to six characters for the names of routines or variables 
has been adopted, using mnemonic definitions. Standard Fortran 77 has been 
adhered to as far as possible. The main exceptions are: 

• The include statement is used for common blocks 

• Input data are read in free format, which is less prone to typing errors 

• Neighbour lists are declared as integer*2 to save memory 

• Continuation lines have '&' in column 6 (non-standard character). 

3.3 Input Parameters 

To be flexible and useful, a general-purpose program should allow a choice of 
input data. Since parameter space is large, the permissible range cannot be 
specified uniquely in advance and every new simulation becomes a test case. 
The complete list of input parameters is defined in Table 2, where each input 
record is separated by a space. The corresponding list of options (denoted 
kz(j)) can be found in routine define. 

Four main types of input may be distinguished. The integration parameters 
etai = 0.03 and etar = 0.06 usually give reasonable results, although slightly 
smaller values are preferable. Further discussions of time-step criteria and 
accuracy in the AC method have been given elsewhere (Makino 1991a, Makino 
& Aarseth 1992). The optimum size of the neighbour sphere depends on the 
problem; usually nnbmax ~ 10 + N 1 ^ 2 is sufficient for small N. For larger N 
(> 1000), Makino & Hut (1988) suggest nnbmax ~ (A^/8) 3/4 . However, see 
Spurzem (1999) for a different point of view. 

An estimate of the initial neighbour radius (rso) for routine nblist should 
be consistent with the scaled density distribution. For homogeneous systems 
(12) gives R s = (2ni/Ny^ 3 Rh, where n\ is the desired membership (say 

2 The code can be downloaded from http://www.ast.cam.ac.uk/~sverre. 
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Table 2 

Input parameters 
Variable Definition 

KSTART Control index (1: new run; >1: restart; > 2: new parameters) 
TCOMP Maximum computing time in minutes (saved in CPU) 



N Total particle number (< nmax) 

nfix Output frequency of data save or binaries (# 3 & 6) 

NRAND Random number sequence skip 

NNBMAX Maximum number of neighbours (< LMAx) 

nrun Run identification index 



etai Time-step parameter for irregular force polynomial 

etar Time-step parameter for regular force polynomial 

RSO Initial neighbour sphere radius 

deltat Output time interval in units of the crossing time 

tcrit Termination time in units of the crossing time 

qe Energy tolerance (restart if de/e > 5*QE and # 2 > 1) 

EPS Softening parameter (square saved in EPS2) 

kz(j) Non-zero options for alternative paths 

xtpari Mass of external Plummer model (# 15 = 1; scaled units) 

XTPAR2 Length scale for Plummer model 

ZMGAS Mass scale for external logarithmic potential (# 15 = 2) 

RGAS Length scale for logarithmic potential 

alphas Power-law index for initial mass function (routine data) 

BODYl Maximum particle mass before scaling 

bodyn Minimum particle mass before scaling 

Q Virial ratio (routine SCALE; q = 0.5 for equilibrium) 

vxrot xr-velocity scaling factor (> for solid-body rotation) 

vzrot z-velocity scaling factor (not used if vxrot = 0) 

rbar Virial radius in pc (for scaling to physical units) 

zmbar Mean mass in solar units (for scaling) 

xcm Displacement of subsystem (routine subsys; # 17) 

ecc Eccentricity of relative motion for subsystem {ecc < 1) 



nnbmax /2) . If necessary, the neighbour sphere is doubled or reduced by the 
volume factor until the membership falls in the acceptable range [1, nnbmax]. 

Another set of parameters specify the initial conditions. The main decision- 
making is controlled by time intervals (deltat, tcrit) and error checking (qe) 
for conservative systems. Among the optional features are rotation, external 
potentials and the creation of two initial subsystems. To guard against typing 
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errors, routine verify checks some input parameters. 



4 NBODY2 Algorithms 

In the following subsections, we discuss some of the main procedures and 
algorithms used by the code nbody2, including several optional features. 

4-1 Initial Conditions 

iV-body simulations usually employ quite specific initial conditions. However, 
for illustrative purposes we include a brief discussion of routine data. 

The present version provides a choice (option 5) between an artificial system 
of constant density and a centrally concentrated Plummer model. Local input 
parameters are first read (exponent a s , maximum and minimum mass m\ 
and mjv)- If ck s = 1 or mi = m^, all the individual masses are taken to 
be unity, otherwise they are selected from a smooth power-law distribution, 
where a s = 2.35 gives a Salpeter-type function. 

In the simplest case (kz(5) = 0), a uniform spherical system is produced by 
choosing x, y, z at random inside a unit sphere with corresponding isotropic 
velocities having randomized magnitudes in [0, 1]. The second choice (kz(5) = 
1) provides for a Plummer model with isotropic velocity distribution (cf. Aarseth, 
Henon & Wielen 1974 for a detailed algorithm). 

4-2 Scaling 

Results of iV-body simulations are often presented in a variety of units which 
hinders comparison with other work. In the present code, a convenient set of 
units are introduced in routine scale; however, all scaling can be bypassed 
if desired (cf. option 16). Once the initial conditions have been generated, all 
coordinates and velocities are first expressed with respect to the centre of mass 
rest frame. 

We introduce the so-called 'standard' units recommended by Heggie & Math- 
ieu (1986), in which G — 1, M = 1, E — —1/4, where M is the total mass 
and Eq is the initial energy. These units are suitable for most types of bound 
systems. This choice corresponds to a virial equilibrium radius R v — 1, rms ve- 
locity V v = y/2/2 and a mean crossing time T cr = 2R V /V V = M 5/2 /(2 |£ |) 3/2 = 
2 a/2. The scaling to internal units is carried out in four separate stages: 
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(1) Scale all the masses by the old total mass, Mq M , to give Mq = 1 

(2) Calculate the total kinetic, potential and virial energies 

(3) Scale the velocities by the specified virial theorem ratio (q) 

(4) Rescale coordinates and velocities to yield the desired total energy. 

Having transformed to scaled units, the output interval (deltat) and termi- 
nation time (tcrit) are re-defined in terms of the crossing time. Using the 
input parameters for the virial radius {rear) in pc and mean mass (zmbar) 
in solar masses, we also introduce scaling factors for conversion to physical 
units. This gives rise to the time unit tstar (in 10 6 yr), velocity unit vstar 
(km s" 1 ) and a re-scaled mass unit zmbar (M ). 

4-3 Integration Algorithm 

We begin by summarizing the starting procedure of the AC method. Given the 
initial conditions r i; Vj and a suitable value of R s , the two force polynomi- 
als are first initialized according to (5) while the neighbour list is constructed. 
We also form the total force (11) and its first derivative, to be used for co- 
ordinate predictions and neighbour corrections. During the second stage, the 
next two orders (6) are calculated for each component, using the neighbour 
list for identification. The initialization procedure is completed by specifying 
irregular and regular time-steps, Ati and ATj, according to (10), whereupon 
the divided differences are formed by (8). The main integration cycle consists 
of the following significant steps: 

(1) Select the next particle, i, to be advanced and define the current epoch 
by t = U + AU 

(2) Make a new sorted time-step list if t > tL and adjust Ati 

(3) Compare t + Ati with To + AT« to decide between a regular force pre- 
diction [case (1)] or a new total force summation [case (2)] 

(4) Predict coordinates of neighbours [case (1)] or all particles [case (2)] 

(5) Combine polynomials for % and predict r t and v t to order F® 

(6) Evaluate the irregular force F? 1 ^ and update the times tk 

(7) Form new irregular differences and include the D 4 term 

(8) [Case (1) only] Extrapolate the regular force and first derivative to give 
F t and F^; proceed to step 15 

(9) Obtain the new irregular and regular force, F™ and F"™, and form a 
temporary neighbour list. Check for n\ = or n\ > n max 

(10) Adjust the neighbour sphere R s and update the times Tfc 

(11) Construct new regular differences and include the D 4 term 

(12) Set F t from (11) and F^ using (4) for both types 

(13) Identify the loss or gain of neighbours and accumulate derivative correc- 
tions; update the neighbour list and convert to differences by (8) 
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(14) Prescribe the new regular time-step ATj 

(15) Update the new irregular time-step At-i 

(16) Repeat the cycle at step 1 until termination or output. 

Here both force polynomials are integrated to the same order; however, two 
separate time-step parameters are used in (10). Note that the decision to re- 
calculate the regular force is based on the next estimated irregular time-step. 
Force polynomials and their derivatives must be combined with care. Thus the 
second regular difference required at step 5 (and at output times) is obtained 
from a differentiation of (4) at t ^ T which yields an extra term in D 3 . 
Likewise for step 8, the regular force and its first derivative are evaluated to 
highest order at an intermediate time before the contributions are added to 
the respective irregular parts. 

Various stratagems may be used to increase the efficiency of determining the 
next particle to be advanced. The earlier scheduling algorithm of searching a 
list of ~ N 1 ' 2 pre-selected members has been modified to include sequential 
sorting. At time t > t L (with t L = initially), a list is formed of all the 
particles due to be advanced during the interval £l + Atj,] and the list is 
then ordered sequentially. 

To include the case of repeated small time-steps and ensure that the global 
time increases monotonically, we have adopted an insert procedure to maintain 
the sequential ordering. Thus the index of a particle satisfying t + Ati < tL, 
evaluated at the end of the cycle, is inserted at the appropriate sequential 
location (routine insert). An updating of the sorted list is made if t > tL, 
whereupon the interval At^ is stabilized on a membership chosen as a com- 
promise between the cost of sorting the quantities tj + Atj and inserting a 
small number of particles in the sequential list. This is achieved by employing 
two variable stabilization factors (in the range 0.25 — 4.0) which are used to 
determine the optimal membership (oc nnbmax). In spite of this complica- 
tion, comparison with the so-called 'heap-sort' algorithm advocated by Press 
(1986) favoured the present algorithm above N ~ 100, although the cost of 
both procedures is oc A r log 2 iV. 

4-4 Force Polynomials 

At several stages, it is necessary to combine force polynomials and their deriva- 
tives. The prediction of the regular force is straightforward, using (2) with 
appropriate time arguments t — Tk (cf. routine nbint). 

The first force derivatives are usually evaluated at the end of an integration 
step of either type, since coordinates or velocities at a general time t ^ to or 
t 7^ T are obtained by an expansion with respect to the interval t — t . In the 
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case of the irregular component, t = to and (4) applies directly; likewise for the 
regular component at an end-point when t = T . However, the prediction of the 
regular force derivative at the end of a sub-step which does not coincide with a 
total force calculation requires special care. Applying one time differentiation 
of (2), we obtain 

Ff } = D 3 (t' t[ + t> t> 2 + i'/ 2 ) + D 2 (t' + t[) + D\ (16) 

where again t' k — to — tk denotes the relevant time intervals (since to = t 
here) . Setting tk = Tk then yields the predicted regular force derivative which 
is added to its irregular counterpart, obtained by (4), to form the total force 
derivative, employed in the subsequent predictions. 

The analogous expression for the second regular force derivative contains an 
extra term t' = to — To in the corresponding coefficient of D 3 (cf. (4)). It is 
required for the prediction of particle % at the beginning of an irregular step 
as well as at output and other high-order predictions. 

Force polynomials based on divided differences result in relatively simple code. 
However, the integration is most conveniently carried out in terms of an equiv- 
alent Taylor series, with (optimized) prediction to order F^ 3 ^ given by 

r t = ((((0.6PV + Y {t) )\-t! + F (1 V + F)t! + v ) t! + r 

v t = (((0.75F (3) t' + F {2) )^' + F (1) )1.5t' + F)2t' + v . (17) 

9 

Here the integration factorials in (9) and factorials in (4) are absorbed in F^ k \ 
and t' — t — to represents the integration interval. 

The coordinate and velocity increments of (17) due to the corrector D 4 contain 
four terms since all the lower derivatives are also modified in (4). Consequently, 
we combine the corresponding time-step factors for r t and v t to yield (in code 
notation) 

s 6 = (((- 1' + s 5 ) 0.6 1! + s 4 ) ^ t! + l - s 3 ) t' 

s 7 = ((0.2t' + 0.25 s 5 )t' + -s A )t' + 0.5 s 3 , (18) 

3 

where again all factorials are absorbed in the force derivatives. The coefficients 
are defined by s 3 = t! x t' 2 t!^ s 4 = t' x t! 2 + t[t' 3 + ^35 s 5 = A + 4 +^j> respectively, 
where the old definition of t' k still applies. The extra factor t' 2 omitted from 
(18) is included in D 4 . The semi- iteration is applied to both polynomials here. 
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However, the regular corrections may not be beneficial if there are large force 
derivatives, as in the case of small softening parameters (i.e. e < R^/N). 

Although the AC scheme is based on the concept of nearest neighbours, there 
are situations when this is less useful. Thus the neighbour sphere determina- 
tion for a distant particle may lead to undesirable oscillations and consequent 
time-step reduction. This problem has been circumvented by including a nom- 
inal small central mass in the irregular force if there are no neighbours inside 
a large neighbour sphere (R s > 50 Rh)] hence adopting rt\ — still permits 
large irregular time-steps. 



4-5 External Potentials 



We have included two optional external potentials, which are treated in a 
similar way. The first case (kz (15) = 1) is the Plummer potential 

GM b 

$1 - ~( r 2 + a 2)l/2> ( 19 ) 



where the background mass Mb and length scale a are specified as input (cf. 
xtpari, xt pari) . The corresponding force is evaluated in routine xtrnli and 
added to the regular component. For initialization, the same routine adds the 
force and first derivative to the regular terms. The appropriate contributions 
are added to the potential and virial energies (routine energy). 

We also include an external logarithmic potential of the type 

2 " (R 9 \og 2 (r/R 9 y (20) 



where M g is the mass scale and R g is the length scale of the continuous mass 
distribution (kz(i5) = 2). A second routine (xtrnl2) is constructed in a similar 
way as above, with C g = M g /R g (denoted cgas). 



4.6 Inelastic Collisions 



An optional procedure (option 12) has been included for the coalescence of 
two particles. The neighbour scheme itself facilitates the determination of a 
collision candidate since a search can be made during the irregular force loop. 
Having identified a close neighbour, the osculating semi-major axis, a, may 
be used in the collision criterion. If accepted, the appropriate global index is 
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saved until the end of the integration cycle. The main steps of the coalescence 
procedure (routine coal) are as follows: 

(1) Improve the collision candidate (j) to order F^ 3 ^ 

(2) Obtain the potential energy with respect to both components 

(3) Form the two-body binding energy E b = —mimj/2a 

(4) Define new mass, coordinates and velocities using the index min(i, j) 

(5) Calculate potential energy with respect to the composite body 

(6) Subtract the two-body energy and tidal potential energy from E t 

(7) Reduce N by 1 and update all common arrays (routine remove) 

(8) Predict coordinates and velocities of neighbours to order F^ 2 ) 

(9) Initialize force polynomials and time-steps for the new body. 

This algorithm maintains fairly good energy conservation with respect to the 
newly updated total energy (Et). Since the common arrays are addressed by 
global indices (rather than pointers), the removal of variables associated with 
a given particle (j) entails a compressing of arrays for all indices k > j, per- 
formed in routine remove. In addition, the neighbour lists require reduction 
by 1 for all index references > j, as well as removal of the specified particle. 



4-7 Escaper Removal 



Routine escape (option 13) contains procedures for the removal of escaping 
particles. A nominal tidal radius, R t , is used as distance criterion, with R t = 
10 R h . The escape algorithm can be summarized as follows: 

(1) Select a candidate (i) for escape if |r;| > 2R t 

(2) Evaluate the binding energy Ei = |v 2 + 0$; accept escaper if Ei > 

(3) Check whether nearest particle is escaping or has large perturbation 

(4) Subtract from the total energy E t and TOj from M t 

(5) Correct the irregular force and first difference for each neighbour 

(6) Reduce N by 1 and update all common arrays (routine remove) 

(7) Re-define coordinates and velocities in the cm. frame and correct E t . 

In a second stage, routine cmcorr (option 18) adjusts all coordinates and 
velocities to the centre-of-mass rest frame. Appropriate modification of the 
total energy involves a correction to the new kinetic energy, and the primary 
integration variables, v and r are updated consistently. 
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4-8 Density Centre 



Models of centrally concentrated systems are most conveniently analysed with 
respect to a well defined cluster centre, in accordance with standard observa- 
tional practice. Following analytical estimates and extensive numerical exper- 
iments, Casertano & Hut (1985) have proposed an operational definition of 
the density centre and corresponding core radius. 

The original definition of the core radius has been modified slightly in order to 
obtain a convergent result using a smaller (central) sample (n ~ N/2). Thus 
the core radius is determined by the rms expression 



where r d denotes the coordinates of the density centre and p; = 15/(47rr|) 
is the density estimator defined with respect to the sixth nearest particle, r 6 . 
We also generalize the density estimator to include the actual mass of the five 
nearest neighbours (Spurzem 1991) in order to be consistent with the original 
expression for equal masses (cf. Eq. (V.3) of Casertano & Hut). The density 
centre and core radius are determined in routine core (option 8) according 
to the schematic algorithm: 

(1) Select all particles inside max (3R C , Rh) as data sample 

(2) Form a list of the six nearest neighbours for each selected particle 

(3) Sort the list of six square distances to find the largest (r@) 

(4) Obtain the total mass (M 5 ) of the five nearest neighbours 

(5) Define individual mass densities pi = 3 M§ / (Ait r\) 

(6) Form the density centre = J2 r i Pi I ' J2 Pi 

(7) Obtain core radius R c and average core density, <p> = J2 Pi/ J2 Pi- 

In the AC scheme, central particles tend to have the largest neighbour den- 
sities, ni/R^; hence an appropriate averaging would also yield a well defined 
density centre. Where relevant, any discussions (and algorithms) involving the 
central distance should be interpreted as referring to the density centre (i.e. 
using |rj — r d \), which is updated at output times. 

4-9 Error Checking 

In standard simulations, energy conservation may be employed as a check to 
safeguard the results from any spurious procedures such as force discontinuities 
or inappropriate program usage. Even systems which include collisions, escape 




(21) 
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and/or analytical external potentials may be treated in this way, since the total 
energy is adjusted accordingly. 

By saving all common variables at specified intervals (usually each output), 
the last interval may be re-calculated with reduced time-steps if the energy 
change exceeds the prescribed tolerance. Since unit 1 (option 1) is reserved 
for common save at termination, unit 2 is employed for this purpose. Thus if 
kz(2) — 1, a common save is performed every output, provided the relative 
energy error de (scaled by the kinetic or potential energy) satisfies de < 5*qe, 
where qe is the specified tolerance. 

Modifications of the accuracy parameters take place if kz(2) = 2. Thus if de 
> 5*qe, the previous common variables are read from unit 2, whereupon 
etai and etar are reduced by a factor of 2, and the current time-steps are 
reduced partially (subject to At; > t — U). A second restart is carried out 
if this procedure does not satisfy the accuracy criterion (using the counter 
ndump). If a restart is not required but de > qe, the integration parameters 
are decreased by (qe/de) 1 / 2 , whereas a small increase is allowed if de < 0.2*qe 
and etai < etao (initial value of etai). 



5 Practical Aspects 

The final section is devoted to practical aspects of A-body simulations, and 
also summarizes some features of a test calculation. 



5. 1 Performance and Accuracy 

The litmus test of a good code lies in its ability to produce satisfactory so- 
lutions in a given time. The CPU time requirement is particularly crucial 
for A-body simulations, which are often designed to make full use of the 
available resources, as in the case of dedicated workstations. In view of the 
time-consuming nature of such calculations, code design must inevitably be 
based on a trade-off between speed and accuracy. 

In the present version, we have made a compromise by defining the basic 
variables r , r t , v , t (and global time) in full double precision, with all other 
variables in standard precision. This gives rise to a faster force calculation 
while the dominant terms of the Taylor series are retained to higher accuracy 
at little extra cost. Assuming a Fortran 77 compiler (or better), the precision of 
the code may readily be changed by suitable declarations, taking care to match 
the array sizes of the saved common blocks (routine my dump). Moreover, on 
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some types of hardware, the CPU time may not depend significantly on the 
choice of precision. 

Actual CPU times are problem dependent and any empirical relation is there- 
fore only a rough guide. In the absence of unique conversion factors for other 
machines, there is also further uncertainty in estimating CPU times. Tim- 
ing tests have been carried out for Plummer models near equilibrium over 
one crossing time with N in the range [100, 2000] and softening parameter 
e = AR V /N. A least square fit including all points gave a relation for the 
computing time per crossing time as T c oc iV L91 with a correlation coefficient 
exceeding 0.9997. Excluding N = 100 increased the exponent to 1.96. This 
result may be compared with an older test for N in [25, 250] which gave a 
value of 1.6 (Aarseth 1985). 

The question of accuracy must also be considered, although it is not an easy 
one. It is well known that the numerical solutions for the point-mass problem 
diverge on a relatively short time-scale (e.g. Miller 1964, Goodman, Heggie & 
Hut 1993). However, in view of the chaotic nature of the iV-body problem, 
it is not justified to aim for the highest possible accuracy at the expense of 
curtailing the integration time or particle number. This is particularly relevant 
when using an approximate model for the interaction potential. 

The most noticeable systematic errors are connected with the integration of 
binaries, and this effect is also present when using a softening parameter. We 
can measure the characteristic error of the two-body motion by studying bi- 
naries in isolation, using the equivalent formulation with one force polynomial 
(program nbodyi). For a binary with eccentricity 0.8 (and no softening) the 
systematic error in semi-major axis per revolution is Sa/a = — 3 x 10~ 5 with 
r] = 0.02 (or -6 x 10~ 6 with rj = 0.01), which corresponds to about 140 (200) 
integration steps for each component. 

5. 2 Softening 

The force law (1) is based on the softened potential —Gmj/(rf- + t 2 ) 1 / 2 which 
was introduced in order to model the interaction between galaxies (Aarseth 
1963). This potential represents the Plummer distribution, where the soften- 
ing size is related to the half- mass radius by ~ 1.3 e. We note that the 
maximum force occurs at r = e/y/2, inside which the interaction tends to a 
harmonic oscillator. 

The introduction of a softened potential is a convenient artifact for modelling 
larger systems by relatively small values of N, thereby suppressing close en- 
counter effects. Considerations based on uniform systems in virial equilibrium 
lead to a close encounter definition r d ~ AR h /N for a 90 degree deflection 
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between two typical particles. On the other hand, the need for sufficient res- 
olution places an upper limit on the amount of softening for a given problem. 
A general discussion of these aspects has been given by several authors (e.g. 
White 1976, Gerhard 1981, Governato, Bhatia & Chincarini 1991). Although 
the present code has been designed for softened potentials, it also works for- 
mally for the case e = as long as there are no critical encounters or persistent 
binaries of high eccentricity. 

In view of the slow asymptotic convergence of the standard softened poten- 
tial, it is desirable to consider expressions with a steeper r-dependence in the 
transition region. One such alternative is 04 = —Gm/(r A + e 4 ) 1//4 which has 
been used for simulations of dwarf galaxies orbiting the Milky Way (Oh, Lin 
& Aarseth 1995). Although the force evaluation is now more expensive, the 
resulting force is a better approximation to Newton's Law for r > e. In the 
AC scheme, the steeper fall-off actually permits the point-mass expression 
to be used for the regular force. Implementation is straightforward, following 
identification of all expressions involving e 2 . 

The interpretation of the softening as an effective half-mass radius has been 
subject to some confusion. Although it has been traditional to generalize the 
interaction potential for two unequal-mass bodies by replacing e 2 with e 2 + e 2 
(e.g. White 1976), the corresponding force does not describe correctly the 
motion of two overlapping Plummer spheres. The actual force law of such 
interactions is rather complicated (cf. Wielen 1979), but alternatives have 
been discussed (Dyer & Ip 1993). However, this still leaves softened potentials 
as a useful tool for suppressing two-body relaxation. For example, collisionless 
simulations of collapsing systems may relate softening to half-mass radius (cf. 
Aarseth, Lin & Papaloizou 1988). 

5.3 Special Features 

It is often desirable to perform large simulations in several stages by saving all 
common variables on disc or tape. This can be done by specifying the run time 
at input (tcomp), together with option 1. Alternatively, on Unix machines the 
command 'touch stop' creates a file (stop) which triggers termination at the 
next timing check (every nmax steps). 

A standard restart is made after reading the saved common file from unit 1 
(using k start = 2). For increased flexibility, some of the input parameters 
can be changed at restart time. Depending on the value of the control index, 
routine modify reads one or two input lines (with the convention that only 
non-zero variables are altered): 

• DELTAT, TNEXT, TCRIT, QE, J, KZ{J) (KSTART = 3 Or 5) 
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• ETAI, ETAR (K ST ART = 4 Or 5). 



Among the various types of initial conditions, we include an example of gen- 
erating two separate subsystems (routine subsys, option 17) which may be 
used for interacting binary experiments. Here the second system is a replica of 
the first (after it has been scaled in the usual way) and the two clumps are as- 
signed initial displacements (±xcm) with specified orbital eccentricity {ecc), 
whereupon the particle number and total mass are updated. In this case, there 
is no well defined density centre, and the modification of the neighbour sphere 
by a radial velocity factor outside the core should be omitted (cf. option 10 
in routine regint) . Likewise, the determination of the density centre and the 
half-mass radius should be modified accordingly (i.e. by using two separate 
procedures), or suppressed altogether. 

The general algorithm of the AC scheme requires current values of all the 
positions at each total force evaluation. However, the frequency of full N co- 
ordinate predictions may be reduced considerably without unduly affecting 
the accuracy (option 14). This is permissible because typical neighbours are 
predicted a number of times at each irregular force summation, thus reducing 
the need for frequent updating of all members. Consequently, the full N pre- 
dictor loop is replaced by a neighbour prediction if the membership ri\ exceeds 
the average neighbour number (updated in kz(u) every output, provided the 
initial value is non-zero). Another reduction of the numerical effort is made 
during initialization by omitting distant contributions to the second and third 
force derivatives in (6) (i.e. R > 3R S ). 

5.4 Data Analysis 

The emphasis of the code is on providing a versatile integration tool, leaving 
the question of data analysis to the individual user. The main facility here is 
the data bank (option 3) which creates a file (unit 3) of the basic variables 
(m, r, v and the identity) for each particle at output times tnext (frequency 
nfix). This data may be read by a separate program, with the first record 
(at, model, nrun) containing the size (at) of the subsequent arrays. 

The general output on unit 6 (routine output) provides a summary of the 
main diagnostics. Most of this output refers to common variables. Among the 
other quantities are the time expressed in crossing times {re) and in 10 6 yr (tg), 
half-mass radius (<r>), maximum neighbour density contrast (cm ax), time- 
weighted neighbour number (<Cn>), centre-of-mass displacements (rcm, vcm) 
and z angular momentum (az). 

Further optional diagnostic output is provided on unit 6 by routine bodies. 
Here an individual line is printed for each particle (using option 9 to control 
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Table 3 

Summary of test calculation 



TIME 


NSTEPI 


NSTEPR 


<NB> 


Q 


<R> 


E 


DE 


0.0 








13 


0.0 


1.87 


-0.25008 


0.000000 


2.8 


44387 


5010 


13 


0.38 


1.27 


-0.25008 


0.000007 


5.7 


227263 


37293 


7 


0.59 


0.76 


-0.25003 


0.000054 


8.5 


353397 


60890 


9 


0.58 


0.66 


-0.25002 


0.000011 


11.3 


476792 


83805 


9 


0.58 


0.65 


-0.25001 


0.000005 


14.1 


601234 


105861 


10 


0.62 


0.61 


-0.25000 


0.000009 



the amount). Likewise, a search for significant binaries is made at each output 
time or at a specified frequency, nfix (kz(6) = 1 or 2). Each output line 
contains the identity of components and their masses, binding energy per unit 
mass, semi-major axis, mean motion, separation, central distance, eccentricity 
and neighbour number. 

Additional information about the radii of specified mass percentiles, including 
the half-mass radius, is provided by routine lagr (option 7) and is printed 
on unit 6 and/or 7, depending on the value of the option. Further analysis 
or plotting routines may be called from routine output, since all the current 
coordinates and velocities (r, v) are known at this stage. 

5.5 Test Calculation 

It is instructive to illustrate the general workings of the code by quoting a test 
calculation, although strict reproducibility may not be achieved on different 
machines. We choose input parameters for a collapsing system of 250 particles 
(cf. test input file). Some selected variables are displayed in Table 3 at every 
other output time, using code notation. 

Comparison of the irregular and regular integration steps (Columns 2 and 3) 
shows a ratio of about 6, for typical average neighbour numbers (Column 4) 
of 10. In spite of the violent relaxation, the total energy (Column 7) is quite 
stable, with fairly small relative energy errors (Column 8). The main error 
occurred at t — 2T cr when the time-step parameters were reduced to 0.012 
and 0.024 (from 0.02 and 0.04), respectively; this was followed by a small 
increase. Note that the half-mass radius (rscale) should be updated fairly 
frequently during phases of violent relaxation (cf. (12)). 

The test calculation required a total CPU time of 35 s on a Sun UltraSparc 1 
workstation, with clock speed 140 MHz and SPEC95 rating of 7.9. Although 
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the AC method gains in efficiency with respect to a single polynomial formu- 
lation as N increases, it can also be used for relatively small systems. Thus 
a second similar example with iV = 25 conserved the total energy to about 
3 x 1(T 5 (using nnbmax = 10 and e = 0.25) after 10389 and 5134 irregular 
and regular steps, respectively. 

5.6 Program Modifications 

The variety of problems suitable for direct iV-body simulations is so large that 
one code cannot deal with all the requirements. However, the basic solution 
method is sufficiently flexible for new procedures to be included. 

Although the present code has been designed for scalar machines, it can readily 
be modified to exploit special features available on vector processors. The 
main speed-up occurs for long loops, where the optimum size is sometimes 
'quantized' to powers of 64 (e.g. the CRAY supercomputer). Thus it may be 
advantageous to increase the size of the neighbour array; further discussion 
can be found elsewhere (Makino 1986, Aarseth & Inagaki 1986). The total 
force loop itself should be simplified by omitting the velocity-dependent part 
which becomes less effective with increasing neighbour number. Even for scalar 
machines, the implementation of a hierarchical time-step scheme would be 
beneficial (McMillan 1986, Makino 1991b). 

Since the square root calculation is by far the most time-consuming part of 
a direct iV-body code, it may be desirable to consider alternative procedures. 
A fast algorithm based on a non-uniform look-up table for 1/r 3 has proved 
effective on scalar machines, where indirect addressing is not a problem (cf. 
Aarseth, Palmer & Lin 1993). This algorithm would also be suitable for the 
regular force calculation without significant loss of precision. 

The code may readily be modified to include a population of massless test 
particles. Now all N particles may be integrated in the standard way, with the 
force loops performed over a smaller number of massive bodies. When mod- 
elling close interactions, it is sometimes desirable to include partially inelastic 
effects. The collision determination would be similar to the coalescence proce- 
dure but the velocity of each component should be modified by the coefficient 
of restitution (Aarseth, Palmer & Lin 1993). 

The effect of an external potential may also be studied as a perturbation of the 
internal motions. This description may be used for star clusters inside a galaxy 
or dwarf galaxies orbiting a large galaxy. In the case of open clusters, the 
assumption of circular motion near the plane of symmetry permits simplified 
expressions for the tidal force (cf. Aarseth 1985). For more general types of 
orbits, it is convenient to employ local rotating coordinates with respect to a 
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guiding centre (Oh, Lin & Aarseth 1992). 

Cosmological modelling by direct N-body methods is limited to much smaller 
particle numbers than can be studied by FFT or tree codes. The standard 
AC code can be used for cosmological simulations as it stands (i.e. q > 1). 
However, it is more efficient to employ comoving equations of motion. Such a 
formulation has been described elsewhere (Aarseth 1985) and is available as 
a separate program called comove. 

Finally, we mention the generation of output for computer movies. This re- 
quires monitoring an additional time-scale after each integration cycle. All rel- 
evant particle coordinates can then be written to a separate file, using scaled 
integer format of two bytes to save memory. The actual movie production de- 
pends on available software; by now it is possible to make nice movies on a 
PC using standard graphics packages. 
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